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Abstract —In this article, the stiffness matrix of a column with linearly variable rectangular cross section along 
the longitudinal axis is obtained. As it is well known in structural analysis, the stiffness matrix can be obtained by 
inversion of the flexibility matrix. So, the terms of this matrix for variable section beam are obtained via the 
energetic method of the principle of virtual works. Parameters called ct^cCf, (If and £ are obtained and they are 
valid only for sections with variable dimensions. The stiffness matrix is then obtained in function of these 
parameters, from de inversion of the flexibility matrix.Finally, the modal analysis of the abutment is carried out 
for the case of the support of a bridge 100 meters high, whose numerical validation of the new bar finite element 
is performed by means of an exact solution by the Continuous Medium Technique (CMT) and by modeled finite 
elements in Ansys academic version. 
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I. INTRODUCTION 

In the structural analysis, the Method of Displacements 
is relevant, which consists of resolving the structure by 
obtaining the initial deformations and, following these 
deformations, it is necessary to draw the diagrams of the 
internal forces of the beam that make up the structure 
(SUSSEKIND, 1978) [1]. The incognitos of the problem 
are: the angle of rotation and the linear displacements at each 
node to be properly constrained from the displacements. 
This should be considered for the various columns. 

As illustrated in Kassimali [2], the advantage of using 
the Method of Displacements over the Method of Forces is 
based on the fact that there is no difficulty in choosing the 
incognitos, since the fundamental problem is only one per 
structure. Thus, the definition of hypergeometric grades 
emerges, which is the number of displacements (linear and 
angular) that can occur at the nodes of a given structure. For 
a high-length Bridge Pillar, it is economically relevant to 
adopt a variable cross section, with a more robust base than 
the top. 

Thus, by adding such a taper in the column shaft it is 
constructively more feasible that such a variation is linear, 
which motivates this article with exact attainment of the 


Stiffness Matrix of a beam element with rectangular linearly 
variable cross section along the axial axis x. In the field of 
approximate solutions in such determination, Fuo et al. [3] 
and Brown [4] are quoted.Thus formulating a New Bar 
Finite Element. 

II. NEW BAR FINITE ELEMENT FOR LINEARLY 
VARIABLE RECTANGULAR SECTION 

a) Method of Forces 

Also called Flexibility Method and developed by James 
Clerk Maxwell in 1864, in order to determine deflections in 
trusses with beam links by labeled joints. To this end, the 
internal energy activated only by the normal effort was used 
[5], Still, based on Timoshenko [5], when promoting the 
temporal delimitation of contributions in structural analysis 
in 1874, Otto Mohr's technique arises in response to the 
problem of nodes (beam joining) of trusses formed by two 
or more rivets,making the connections no longer labeled 
and, therefore, the need to compute the deformation internal 
force energy with all the internal forces. In equation 1 (a), 
the calculation of deflections is verified by the Maxwell 
formulation for labeled trusses. On the other hand, in 
equation 1 (b), the same occurs, however, for connections 


www.iiaers.com 


Page | 59 







International Journal of Advanced Engineering Research and Science (IJAERS) 
https://dx. doi. orq/10.22161/ijaers. 75.8 


[Vol-7, Issue-5, May- 2020] 
ISSN: 2349-6495(P) / 2456-1908(0) 


with two or more fixations per node, according to Mohr's 
formulation. 
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with: S i0 = Nodal deflexion to the t node; N i: M b V L — Axial 
load. Flexion moment. Bending moment and Shear force for 
the real load of the structure; and N 0 , M 0 , V 0 = Axial load, 
Bending moment and Shear force for the virtual load via 
P.V.W. acting on Node i. 

Yet, in Charlton [6], Maxwell's publication entitled: “On 
reciprocal figures, frames and diagrams of forces” it is 
evident. In the work of the same author, page 83, Mohr's 
publication entitled: " Beitrag Zur Theorie Du 
Bogenfachwerk trager " can be seen. In 1886, Heinrich 
Muller - Breslau postulated the systematization of the 
Method of Forces defined earlier by Maxwell and Mohr 
(KINNEY, 1957) [7], The publication of such 
systematization is: "Die neue methoden der 
baukonstruktionen". The fundamental basis of the Method 
of Forces is the compatibilization of the angular and linear 
displacements of the extracted connections in order to 
establish extra equations and make the system of equations 
of the problem linearly independent and, therefore, make it 
solvable. 

The compatibilization of the displacements due to the 
connections extracted via definition of said method will be 
adopted from the positive convention for linear 
displacements in the positive direction of the jr and y axes 
and rotation with vector notation in the same direction as the 
positive z axis. Thus, using a formulation present in Kiseliov 
[8], the system of displacement compatibility equations is 
written in matrix notation as: 

{d} = [F].{Y} + {S} + {5 r } (2) 

with: {d} — Vector of displacements in the initial 
hyperstatic structure; {5} = Vector of displacements in the 
fundamental problem; {S 1 } = Vector of displacements in 
the thermal problem; {X} — Vector of incognitos for the 
Method of Forces; and [F] = Flexibility Matrix. 

In order to assemble the Flexibility Matrix [F]. just 
consider the derived systems to become, in the background, 
both the fundamental problem responsible for the analysis 
of the actual and active load in the initial structure and the 


thermal problem. In addition, the vector {5} is neglected in 
this analysis, since it is only desired to obtain the Flexibility 
Matrix. Such Systems derived from a fixed-fixed beam are 
shown in Figure 1 as much as the Main System. 


Mi = 1 M f = 1 j 

lT\ *r = 1 (T\ 



Fig. 1: Method of Forces: (a) Main System, (b) 1st. 
Derived System, (c) 2nd. Derived System, (d) 3rd. Derived 
System 

At the initial node there is the Degree of Freedom 0 t . The 
final node is characterized by the degree of freedom also in 
rotation Of and, finally, Sf in linear displacement for the 
final node. The system of equations is expressed by: 
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with: cq, a.f — Angular displacement for the initial and final 
node of the beam, respectively; £ — Angular displacement 
in a node contrary to where a t and ctf occur; M t , Mj = 
Flexion moment along the beam with unitary moment 
imposition; V L , Vj- = Shear force along the beam with unitary 
and vertical load imposition; and Nf — Axial load of the 
beam with unitary and horizontal load imposition, 
b) Obtaining the displacements a, (3 and £ via PVW 

The Principle of Virtual Works (PVW) was postulated 
by John Bernoulli in 1717 and is based on the Principle of 
Energy Conservation. This principle was also linked to the 
concept of virtual displacement. In virtual displacement, 
when the material point is in equilibrium, real displacement 
cannot occur. And for a particle to be in equilibrium, the 
condition of nullity to the work of all external forces must 
be satisfied (STAMATO, 1983) [9]. Figure 2 shows the 
linear and angular displacements (f, andcq, Uf and £, 
respectively, as: 
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Fig. 2: Internal Force Diagrams: (a)Mi = (b) Mf = 

M f , (c) V = V, (d) N f = N f 

For the beam element of length L and dimensions in the 
cross section H y (x) and Hjx), as shown in Figure 3, the 
cross section area A (X) and the Moment of Inertia 40) 
around the z axis, is written as: 



b z = H z (x = 0) h z = Hjx = L) 

by = Hy( X = 0) ky = Hy(.X = L) 


Fig. 3: Beam Element with dimension in the linearly 
variable cross section along the axial axis x 

40) = k 4 . x 4 + k 2 . x 3 + k 3 . x 2 + k 4 .x + k 5 (5 a) 

A(x) — k 6 .x 2 + k 7 .x + k 8 (5 b) 

where: k 1 = AC 3 ; k 2 = C 2 .(3.A.D + B.C ); k 6 = A.C ; 
k 3 = 3. C. D. ( A. D + B. C); O = B. D 3 ; k 8 = B. D; 
k 4 = D 2 . (A D + 3. B. C. £>); k 7 = A. D + B. C; 


B — b z ] D — b v ] A — 


h 7 — b y 


C = 


with: b z , h z —Cross section dimensions, parallel to axis z, in 
the initial and final section, respectively; b y , h y —Cross 
section dimensions, parallel to axis y, in the initial and final 
section, respectively; E —Longitudinal Elasticity Module; 
G — Cross Elasticity Module; 40 ) — Variation along the 
axial axis x of the inertia moment surrounding axis z; 
A(x) — Variation along axial axis x of the cross section 
area; k c — Shape factor and L — Beam length. 


Analysing the Area Static Moment Q(x) for a 
rectangular section with dimensions H y O) and H z (x), 
characterized in Figure 4, results in: 
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2. H y (x) 

l(x) — A'(x).z'(x) — 

Proceeding the calculation of the shape factor k c , after 
insforming the integration in area A into along the length, 
ere is: 

A Q) J Q 2 O) 
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( 7 ) 


When applying equations 6 (a - d) in equation (7) and 
ying out the integration in z and consequent 
plifications, we conclude that the shape factor A: c remains 
hanged along axis x with the following value: 

5 

= 6 (8) 
Finally, when applying the equations (5), 6 (a - d) and 
in equations 4 (a - d), the flexibility coefficients cq, cry-, 
id (3f are attained, expressed by: 

> if A. D A B.C. 
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where A, B, C and D are parameters linked to geometrical 
dimensions of the cross section. And yet, list yourself: y = 
L.(A.D + B. C). (2. B. D + 3 .A.D.L - B.C.L); 

T] 3 = -2. D 2 . L 2 . (A. D +B. C) 3 ; J? 4 = l 2 . (A. D + B. C); 

2. % = -L. (A. D-B. C). (2. B.D +A.D.L + B. C. L); 

77s = -L. (A. D-B. C). (2.B.D -A.D.L + 3. B. C. L); 

? 7 2 = 2. D 2 . (A. L + B) 2 ; ? 7 6 = 2. B 2 . (C. L + D) 2 ; 

77-7 = 2. L 2 . (C. L + D) 2 . (A. D-B. C) 3 ; 

? 7 10 = D. L 2 . (C. L + D). (A. D + B. C) 3 ; 

? 7 9 = B.D.(A.L + B).(C.L + D). 


c) Method of Displacements 

Such method consists of resolving the structure by 
initially obtaining the deformations by means of the internal 
forces on the beams that make up the structure, while in the 
Method of Forces, extra bonds are extracted to structural 
staticity and stability. In the Method of Displacements, the 


locking of the bound nodes of the structure is promoted in 
order to attain the connection of the beams by means of 
embedments. In matrix notation, the system of equations of 
equilibria of the unbalancing forces by nodes is expressed as 
follows: 


{M} = [K].{D} + [y] + {y T } (10) 

where: {M} = Vector of unbalancing forces; [K] = Stiffness 
Matrix; {D} = Vector of incognitos in the Method of 
Displacements; {y} = Vetor of unbalancing forces in the 
fundamental problem fundamental; {y 7 } = Vector of 
unbalancing forces in the thermal problem; and 
fcj, kf, a, Tf — Stiffness coefficients. 


Using the definitions in Kassimali[2] and Kiseliov [8], 
the Stiffness Matrix [K] will be the inverse of the Flexibility 
Matrix [F]. The terms of stiffness are concluded, such as: 
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The present formulation was proposed by George Alfred 
Maney in 1915, and was also denominated Rotation-Arrow 
Method. See [10]; [11] and [12], 


HI. MATRIX CONDENSATION 

Based on Paz [13] and the application in the calculation 
of natural frequencies in Alves Filho (p. 200) [14], matrix 
condensation consists of rewriting the system of equations 
in terms of some of its variables. In Figure 5 shows a pillar 
with n subdivisions and (n + 1) nodes, as well as the 
resulting degrees of freedom (S and 9) and nodal forces, in 
order to exemplify such a condensation procedure. 
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Fig. 5: column subdivision by finite elements of bar 

The ODE’s system remaining, for the non-damped 
vibration, expressed by: 

0)1 


[Mqq] [ Mgs] 
[M S g] [M SS ] 
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{F} T = {F 1 F 2 F 3 ... F r 

The submatrix [Mg g] presents the terms of rotational 
masses, with little representativeness (in magnitude) in 
relation to translational masses (included in the submatrix 
[Mgs]). From this statement, sub-matricesfMg^] and [Msg] 
are also disregarded, and equation is rewritten (12) as: 


[ 0 ] [ 0 ] 

[o] [Mss] 


W 

' It*} 


+ 


[Kgg] IXesl 
[Kse] [Fssl 


f{0}) _ f{M}) 

1{F}J 


(13) 

Expressing equation (13) in the form of a matrix 
equation system, we have: 

[Kgg].{8} + [Kg S ].{8} = {M} (14 a) 

[M ss ]. {5} + [K se ]. { 6 } + [K ss ]. {5} = {F} (14 b) 

In order to express the displacements { S }, the rotation 
vector { 9 } is isolated in equation (14 a) and applied to eq. 
(14 b), concluding: 

[Mssl{8} + [K S g].([Kgg]-\{M] - [Kgg]-\ [/f efi ].{<5}) 

+ [KssU8} = {F] (15) 

rearranging, there is: 

[M ss ].{ 5 } + [/r].{ 5 } = {F*] 

With: [F‘] = [K SS ] - ([K S gl [Kgg]~\ [Kg S ]f, 

{F*} = {F} 


(15 a) 


v sst 

(.[KsgliKggrfiM]). 


and: \I(*] the condensed stiffness matrix, and {F*} the 
condensed vector of transverse forces. 

Condensed ODE for proportionally dampened 
vibration([C*] = a m . [Mss] + a k- [F*])is expressed by: 

[M SS l {5} + [C*]. {8} + [K*]. {5} = (F*} (16) 

IY. VERIFICATION OF (FEM) MODELING VIA 

(CMT) 


Through the dynamic analysis of the pillar by the 
Continuous Medium Technique (CMT) processed in Melo 
[15], the differential equation of the problem is written, as: 

-[/]■ {v""} + [5], {v "} + [M].{v] = {V/} (17) 

Based on the WallPanels Theory (WPT), the differential 
equation of the dynamic request of the massive column 
shown in Figure 6 is expressed, through equation (17), as: 

~[J]-{q""(.x,t)} + [S].{q"(x,t)) + [M].{q(x, t)} 

= {V/(x,t)} (18) 

with: [/] is the column stiffness matrix; [5] is the lintel 
stiffness matrix (and for the massive pillar, modeled in this 
item of the thesis, it will be null); [M] is the mass matrix of 
the abutment and q ( x, t) is the function of the displacements 
dependent on space and time. Through the harmonic 
analysis of equation (18) and imposing the procedure for the 
separation of variables, it is written: 

—j.u""(x) —m.g(t ) 


= ~An 


(19) 


u(x) g(t ) 

with: q(x,t ) = u(x).g(t). The characteristic equation of 
ODE written in space is expressed, via equation (19), as: 


j. a> 4 - X a = 0 

and by equation (20) solution is expressed: 


Cl) = 




( 20 ) 


( 21 ) 


By calculating the stiffness, where j = E.I, together with 
Pfeil (p. 211) [16], concludes by the rigidity of the column 
shown in Figure 6, the following: 

j=E.I = E.fU Ztopo 

, 10 m . (5 m) 3 

= 4,67 x 10 4 Pa . 1,2793 .---— 

12 

„N 

= 6,22326 x 10 12 — 

772 

( 22 ) 

As presented inDziewolski[17]followed by the 
adjustment of the column stiffness using the coefficienta = 
Vl 5 (f° r si m pl e structures)and is expressed via equation: 


Jf= 


H 2 


. a 


6,22326 X 10 12 — i 

= 105,29633 


(100 m) 2 ' 1,5 


(23) 
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Applying to equation (23) in equation (21) the first 
vibration frequencies of the bridge pillar, shown in Figure 6, 
are written as: 


105,29633 


(24) 


Y. COLUMN MODELING IN 5 FE APPLIED TO 
MODAL BRIDGE ANALYSIS 

In order to exemplify the use of the linearly variable 
rectangular section stiffness matrix in bridge pillars (See 
Figure6), cross-sectional dimensions at the base b y = 
12.5 m and b z = 25 mand, at the top, h y = 5 m and h z = 
10 m, and the modes of vibration are obtained through 
modal analysis [18], modeling via five finite elements [19] 
and general formulation for n mass presented in Warburton 
[20] and matrix condensation [21]; and [22]. The material 
used in the bridge is reinforced concrete of resistance class 
C - 40 [23]. Therefore, the Longitudinal Elasticity Module 
will be E — 35 xlO 9 Pa and the Poisson Coefficient will be 
v = 0,20.In order to validate the example, modeling is 
performed using the ANSYS academic version software. It 
should be noted that for the vibration modes, the dimensions 
adopted for the cross section become irrelevant. 



Fig. 6: Bridge with linearly variable section columns 

From this analysis, the coefficients for the generation of 
the column stiffness matrix and vibration frequencies uij and 
autoversors Aj, via the nullity of the determinant | [ K\ — 
A. [M] |, are obtained in Table 1, as follows: 

A\ = 6.16906 x 10 s (f~~j : A. 

_ / rad\ 2 

= 7.14143 x 10 s ( —J ; 


32 _ 


_ irad\~ 
7.83961 x 10 s (- ; 


V 5 ) 


/rad\ 


= 8.34502 x 10 s ; 


Al = 8.71653 x 10 s 


A 2 i 


(?) 2 


rad 


= 7.85434 -= 1.25006 Hz; 

H z s 


rad 

s 

rad 


o) 2 = 8.45070 -; a> 3 = 8.85415 

s s 


rad 

s 

rad 


(25. a — e) 


cl> 4 — 9.13511-; cn 5 = 9.33624 - (26.a-e) 


Table. 1: Parameters for the generation of the stiffness 
matrix of the linearly variable section column 


Bar Finite Element 

1 

2 

3 

4 

5 

A (ADM) 

-0.015 

-0015 

-0.015 

-0.015 

0.015 

B (m) 

25.000 22.000 19.000 16.000 13.000 

C (ADM) 

- 0.075 

0.075 

- 0.075 

0.075 

0.075 

D(m) 

12.500 11.000 

9.500 

8.000 

6.500 

= a f (x 10 s ) 

6.540 

4.430 

2.890 

1.799 

1.050 

k t (x 10 13 ) [N.m] 

2.274 

1.364 

0.759 

0.382 

0.166 

k f (x 10 13 ) [N.m] 

2.261 

1.356 

0.754 

0.379 

0.165 

a (x 10 13 ) [N.m] 

1.134 

0.680 

0.378 

0.190 

0.083 


When considering the five finite elements, the Stiffness 


[If] and Mass [M] Matrices for the pillar are worth: 


[K] 



1.701 x 10 13 

7.634 x 10 11 

3.178x10 

L0 1.196 xlO 9 

3.891 x 10 7 


7.634 x 10 11 

1.019 xlO 13 

4.241x10 

L1 1.596 xlO 10 

5.193 x10 s 

= 

3.178 x 10 10 

4.241 x 10 11 

5.666 x 10 

12 2.132 xlO 11 

6.939 x 10 9 


1.196 x 10 9 

1.596 xlO 10 

2.132x10 

11 2.848 xlO 12 

9.270 x 10 10 


■ 3.891 x 10 7 

5.193 x 10 s 

6.939x 10 

9 9.270 xlO 10 

1.238 x 10 12 - 

[M] 







1.955 xlO 7 

9.995 xlO 5 

0 

0 

0 



9.995 xlO 5 

1.227 xlO 7 

6.761 x 10 5 

0 

0 


= 

0 

6.761 xlO 5 

7.281 x 10 6 

4.374 x 10 5 

0 

kg 


0 

0 

4.374 x 10 5 

4.030 xlO 6 2.685 x 10 5 



0 

0 

0 

2.685 x10 s 1.984 x 10 6 . 



a) Validation via ANSYS 


Figure 7 shows a group of vibration modes, for the 
massive column of the bridge shown in Figure 6, modeling 
in the ANSYS software. 62,468 nodes and 13,635 finite 
elements were used, producing a mesh with 93.99%. 

The validation of the first vibration frequency f t (using 
equation 26 a) is processed and results in: 

7.85434 — 

fi=rr = -^-— = 1.25005 Hz (27) 

An An 
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Comparing the first vibration frequency by modeling in 
ANSYS, see Figure 7, with the value presented in equation 
(27), an approximation of: 


(1.25005 - 1.23300) Hz 
1.25005 Hz 


.100% = 1.36% 


(28) 


B: Pilar Variavel - C-90 

Total Deformation 
Type: Total Deformation 
Frequency: 1,233 Hz 
Unit: m 

12/12/201922:46 


_ 0.00048451 Max 

■ 0,00043068 

— 0,00037684 

- 0,00032301 

- 0.00026917 

— 0,00021534 

— 0,0001615 

- 0,00010767 
_ 5,3835e-5 

■ OMin 


.3 




9000 (m) 


(a) 


B: Pilar Variavel - C-90 

Total Deformation 3 
Type: Total Deformation 
Frequency: 4,5579 Hz 
Unit m 

12/12/2019 22:48 

0,00057471 Max 

0,00051085 
0,00044699 
0,00038314 
0,00031928 
0,00025542 
0,00019157 
0,00012771 
6,3856e-5 
OMin 





B: Pilar Variavel C-90 

Total Deformation 6 
Type: Total Deformation 
Frequency: 10,445 Hz 
Unit m 

12/12/201922:49 


■ 0,00059576 Max 

0,00052956 
— 0,00046337 
- 0,00039717 
— 0,00033098 
— 0,00026478 
- 0,00019859 
- 0,00013239 


6,6195e-5 

OMin 



(b) 


(c) 


Fig. 7: Vibration modes in bridges with linearly variable 
section columns: (a) I s ' mode, (b) 2 nd mode, (c) 3 rd mode 


The 1.36% discrepancy between the fust vibration 
frequency, via manual calculation by finite bar elements and 
by modeling in the ANSYS software, is due to the small 
number of bar elements used in manual modeling. However, 
the formulation present here is quite satisfactory to verify 
the order of magnitude of the results obtained via modeling 
in commercial software. 


b) Validation via CMT 


Using the first root of the polynomial presented in 
equation (25a), the first vibration frequency of the abutment 
(via CMT) is expressed as: 


V6.16906 x 10 5 rad/s 
105.29633 


rad 

7.45927 - 

s 


(29 a) 


and in fundamental frequency, there is: 

a), 

ft == 1.18718 Hz (29 b) 

2n 

concluding by divergence in relation to the modeling 
performed in ANSYS, see Figure 7, the following: 


I frMC ~ fANSYsl 
frMC 


. 100 % 


|1.18718 — 1.2331 Hz 
1.18718 Hz 


. 100% = 3.87 % 


(29 c) 


It is observed that the percentage difference by the CMT 
was greater than by the bar FEM, this due to the 
interpolation performed in the coefficient/? of equation (22). 
As well as, it is verified the use of the adjustment coefficient 
a in equation (23). While in equation (28) the analysis is 
processed by the finite element method, eliminating the 
imposition of such a coefficienta. Even so, it is possible to 
satisfactorily validate the modal analysis of the column 
shown in Figure 6, both by FEM and CMT. 


VI. CONCLUSION 

In this article, the terms of the Flexibility Matrix for a 
variable section were obtained through the internal energy 
activated by the internal forces acting on the beam element 
for the Derivative Systems. The contribution of this 
publication is the obtention of the stiffness matrix by 
parameters easily obtained from the flexibility matrix which 
takes into account the variation of the dimensions of the 
element cross section. When analyzing the shape factor k c , 
it was found to be a constant value according to the section 
shape. After determining the terms of the Flexibility Matrix, 
the terms of the stiffness matrix are obtained as shown in 
this work. 

Finally, it is made explicit that the expressions presented 
here for cq, a.f, £ and (if are valid only for cases of variable 
cross-section along the longitudinal axis of the pillar. In case 
of a constant section, one should resort to the expressions 
widely postulated in technical texts, under penalty of 
mathematical indetermination by these explicit expressions 
for a, (3, £ and [if. 

This contribution emphasizes the achievement of the 
exact stiffness matrix for such a cross-section configuration, 
and the application of such a matrix is performed in the 
calculation of the vibration frequencies of the pillars of a 
bridge with a tray positioned at 100 meters, this relative to 
the base, of the mentioned columns. 

For the first vibration frequency, with the column being 
subdivided into five finite bar elements, an approximation 
of 1.36% is concluded. Such an approximation is excellent. 
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due to the number of finite elements used in the 

discretization of the lOOmeters high column, the validation 

resulting from modeling in ANSYS. 
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